Simulation of cohesive head-on collisions of thermally activated nanoclusters 
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Impact phenomena of nanoclusters subject to thermal fluctuations are numerically investigated. 
From the molecular dynamics simulation for colliding two identical clusters, it is found that the 
restitution coefficient for head-on collisions has a peak at a colliding speed due to the competition 
between the cohesive interaction and the repulsive interaction of colliding clusters. Some aspects 
of the collisions can be understood by the theory by Brilliantov et al. (Phys. Rev. E 76, 051302 
(2007)), but many new aspects are found from the simulation. In particular, we find that there 
are some anomalous rebounds in which the restitution coefficient is larger than unity. The phase 
diagrams of rebound processes against impact speed and the cohesive parameter can be understood 
by a simple phenomenology. 
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INTRODUCTION 

Inelastic collisions are the process that a part of ini- 
tial macroscopic energy of colliding bodies is distributed 
into the microscopic degrees of freedom. This irreversible 
process of head-on collisions may be characterized by the 
restitution coefficient which is the ratio of the normal 
rebound speed to the normal impact speed. Although 
it was generally believed that the restitution coefficient 
is a material constant, modern experiments and simu- 
lations have revealed that the restitution coefficient de- 
creases with the increase of impact velocity. 0, 0, H, HI For 
example, in the case of collisions between icy particles, 
we can easily find the monotonic decrease of the restitu- 
tion coefficient against impact velocity without any flat 
region. The dependence of the restitution coefficient 
on the low impact velocity is theoretically treated by the 
quasistatic theory 0, 0, @, @, E3]. In Ref.jl], Kuwabara 
and Kono performed impact experiments by the use of 
a pendulum of various materials to validate their theo- 
retical prediction. On the other hand, the dependence 
of the restitution coefficient on the high impact veloc- 
ity is treated by the dimensional analysis based on plas- 
tic collisions. ll| From the dimensional analysis, the re- 
lationship between the restitution coefficient e and the 
impact velocity V becomes e oc V^ 1 ^ 4 , which coincides 
with experimental results by the use of a steel ball and 
blocks of various materials such as hard bronze and brass. 
P, d, [ll[ We also recognize that the restitution coefficient 
can be less than unity for head-on collisions without any 
introduction of explicit dissipation, because the macro- 
scopic inelasticities originate in the transfer of the energy 
from the translational mode to the internal modes such 
as vibrations. 0, 0, 

Although it is believed that the restitution coefficient 
for head-on collisions is smaller than unity, the resti- 
tution coefficient can be larger than unity in oblique 
collisions. 0EE (3 For example, Louge and Adams ob- 
served such an anomalous impact in which the restitution 
coefficient is larger than unity in oblique collisions of a 



hard aluminum oxide sphere onto a thick elastoplastic 
polycarbonate plate in which the restitution coefficient 
increases monotonically with the increase of the magni- 
tude of the tangent of the angle of incidence. [3] They 
explained that this phenomena can be attributed to the 
change in rebound angle resulting from the local defor- 
mation of the contact area between the sphere and the 
plate, which causes the increase in the normal compo- 
nent of the rebound velocity against the collision plane. 
The present authors performed a two-dimensional im- 
pact simulation with an elastic disc and an elastic wall 
consisted of nonlinear spring network to reproduce the 
anomalous impacts. They also explained the mechanism 
to appear large restitution coefficient based on a simple 
phenomenology by taking into account the local surface 
deformation. [15j] 

The static interaction between macroscopic granular 
particles is characterized by the Hertzian theory[l7|, G3] 
of the elastic repulsive force as well as the dissipative 
force which is proportional to relative speed of colliding 
two particles. The total force acting between granular 
particles in contact is assumed to be a combination of the 
elastic repulsive force and the dissipative force in the qua- 
sistatic theory, with which many aspects of the inelastic 
collisions for such granular particles can be understood. 
This theory can reproduce the restitution coefficient as a 
function of the colliding speed observed in experiments 
and simulations. @, [H, [3] 

Although the repulsive interaction becomes dominant 
for collisions of large bodies, cohesive interactions such as 
van der Waals force and electrostatic force play impor- 
tant roles for small clusters of the nanoscale. 20, 21. 22\ 
Recently, Brilliantov et al. have developed the quasistatic 
theory for inelastic collisions to explain the relation be- 
tween the colliding speed and the restitution coefficient 
for cohesive collisions. [23| The result of an experimen- 
tal result of collisions of macroscopic particles with the 
cohesive interaction is consistent with the theory. 24 1 

For molecular dynamics simulations of small clusters, 
many empirical potentials are used to mimic the interac- 
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tion between various atoms. [25j Among them, most com- 
monly used one is the Lennard-Jones potential: 



(1) 



which well approximates the interaction between inert 
gas atoms such as argons. [2^, [26| Here, r*y is the distance 
between two atoms labeled by i and j, respectively, e and 
a are the energy constant and the characteristic diame- 
ter, respectively. In this potential, the second term on 
the right hand side represents the cohesive interaction 
which is originated from van-der Waals interaction. 

Dynamics of nanoclusters are extensively investigated 
from both scientific and technological interests. There 
are a lot of studies on cluster-cluster and cluster- 
surface collisions based on the molecular dynamics 
simulation. 27, 28, 2^, 3^, 3l| We observe variety of re- 
bound processes in such systems caused by the compe- 
tition between the attractive interaction and the repul- 
sive interaction of two colliding bodies. Binary colli- 
sions of identical clusters cause coalescence, scattering, 
and fragmentation depending on the cluster size and the 
impact energy. [53, [28[ On the other hand, cluster- 
surface collisions induce soft landing, embedding, and 
fragmentation. [31| The attractive interaction plays cru- 
cially important roles in such colliding processes. 

However, the attractive interaction may be reduced 
in the case of some combinations of the two interact- 
ing objects and the relative configuration of colliding 
molecules [32]. Awasthi et al. carried out the molecular 
dynamics simulation for collisions of Lennard-Jones clus- 
ters onto surfaces to simulate the collision of a Bi cluster 
onto a SiC>2 surface. [33[ They introduced a cohesive pa- 
rameter to characterize the magnitude of attraction and 
investigate the rebound behavior of the clusters. Simi- 
larly, recent papers have reported that surface-passivated 
Si nanoclusters exhibit elastic rebounds on Si surface due 
to the reduction of the attractive interaction between the 
surfaces. [U HH These results suggest that nearly repul- 
sive collisions really exist even in small systems. 

In the case of purely repulsive collisions between two 
identical nanoclusters, we have already reported that the 
relation between colliding speed and the restitution co- 
efficient may be described by the quasi-static theory for 
inelastic impacts, though the restitution coefficient ex- 
ceeds unity for small impact speed. H, 37 1 In addition, 
on the basis of the distribution function of macroscopic 
energy loss during collision, we have shown that our nu- 
merical results can be approximated by the fluctuation 
relation for inelastic impacts. [3^ 

The aim of the present paper is to study statistical 
properties in binary head-on collisions of identical nan- 
oclusters. In particular, we numerically investigate the 
effects of attractive interaction on the restitution coef- 
ficient in rebound processes. The organization of this 



paper is as follows. In the next section, we introduce our 
numerical model of colliding nanoclusters and the setup 
of our simulation. In Section III, we summarize the re- 
sults of our simulation. In Section IV, we mainly discuss 
the system size dependence of our results. In Section V, 
we summarize our results. Appendices A, B, and C treat 
the calculation of the surface tension, the technical cal- 
culation on the system size dependence of the restitution 
coefficient, and stability of spherical shape of a elastic 
droplet, respectively. 



MODEL 

Let us introduce our numerical model. Our model con- 
sists of two identical clusters, each of which is spheri- 
cally cut from a face-centered cubic (SC-FCC) lattice of 
"atoms" . We typically use 682 atoms systems which are 
13 layers SC-FCCs. The system size dependence will be 
discussed in Section IV. Here, we list the relation be- 
tween the number of "atoms" and the number of lay- 
ers in one cluster in Table HI The clusters have facets 
due to the small number of "atoms" (Fig. []}. All 
the "atoms" in each cluster are bound together by the 
Lennard-Jones potential U(rij) in Eq.JT]). When we re- 
gard the "atom" as argon, the values of the constants be- 
come e = 1.65 x 10 _21 J and a = 3.4A, respectively. [25| 








FIG. 1: (Color online) A typical situation of our simulation 
of two colliding clusters. Each of them contains 682 "atoms" 
which are bound by the Lennard-Jones potential. 

Henceforth, we label the upper and the lower clusters 
as C u and C , respectively. We assume that the interac- 
tive potential between the atom k on the lower surface 
of C u and the atom I on the upper surface of C l is given 
by 



(2) 



where r k i is the distance between the surface atom k 
and atom I. We introduce the cohesive parameter c to 
characterize the attraction between the atoms of different 
clusters. 1 33] 
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TABLE I: Relation between numbers of layers and atoms. 



Number of Layers 



Number of Atoms 



3 
5 
7 
9 
11 
13 
15 
17 



12 

42 

135 

236 

433 

682 

1055 

1466 



The procedure of our simulation is as follows. As the 
initial condition of simulation, the centers of mass of C u 
and C l are placed along the z-axis with the separation 
a c between the surfaces of C u and C l . The initial veloc- 
ities of the "atoms" in both C u and C l obey Maxwell- 
Boltzmann distribution with the initial temperature T. 
The initial temperature is set to be T = 0.02e in our 
simulations. Sample average is taken over different sets 
of initial velocities governed by the Maxwell-Boltzmann 
velocity distribution for "atoms" . 

To equilibrate the clusters, we adopt the velocity scal- 
ing method [H, [39} and perform 2000 steps simulation 
for the relaxation to a local equilibrium state. Here let 
us check the equilibration of the total energy in the ini- 
tial relaxation process. Figure EJa) is the time evolution 
of the kinetic temperature of C u , where T denotes the 
scaled temperature by the unit e. This figure shows the 
convergence of temperature to the desired temperature 
T = 0.02e. On the other hand, Figl^b) is the probability 
density distribution of speed of "atoms" in C u when the 
equilibration process is over, where v denotes the scaled 
velocity for "atom" by the unit \J e/m. The solid curve 
in Fig(2]Jb) shows the probability density distribution of 
speed v% of "atoms" indexed by i in equilibrium, 



\{vi) = 4tt 



(, 



m \ 3 / 2 2 

) v l exp 



V27rfcT 



with T = 0.02e. This agreement shows that the upper 
cluster C u is equilibrated during the equilibration pro- 
cess. 

After the equilibration, we give translational velocities 
to C u and C l to make them collide against each other. 
The relative speed of impact ranges from V = 0.02^/e/ m 
to V = 0.6^/e/m. Here, the characteristic speed is the 
thermal velocity for one "atom" y/T/m, where m is the 
mass of each "atom". This situation might correspond 
to the sputtering process or collisions of interstellar dusts 
or atmospheric dusts. Although it is not easy to control 
the velocity of colliding clusters in real nanoscale exper- 
iments, the effects of thermal fluctuation to the center 
of mass of each cluster might be negligible if clusters are 
flying in vacuum. 
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FIG. 2: (a) Relaxation of kinetic temperature and (b) distri- 
bution of speed of atoms after equilibration process. T and 
v are temperature and speed of "atom" scaled by units e and 
y/e/m, respectively. 



Numerical integration of the equation of motion for 
each atom is carried out by the second order symplectic 
integrator with the time step dt — 1.0 x 10~ 2 <7 / \Je/m. 
To reduce computational costs, we introduce the cut-off 
length a c = 2.5a of the Lennard- Jones interaction, which 
sometimes affects the energy conservation of a system 
although the Hamiltonian of the system is conserved. We 
have checked that the rate of energy conservation, \E(t) — 
Eq\/\E \, is kept within 10 -5 with the cutoff length a c = 
2.5cr, where E is the initial energy of the system and 
E(t) is the energy at time t. In general, the value between 
3er < <j c < 4ct is used for the energy conservation about 
\E(t)-E Q \/\E \~lQ-*. 

We let the angle around z— axis, 9 Z , be 9 Z = when 
the two clusters are located mirror-symmetrically with 
respect to z — 0. In most of our simulation, we set 9 Z 
at 9 Z — as the initial condition. The dependency on 9 Z 
will be shown in the next section. 



RESULTS OF OUR SIMULATION 



Relation between impact speed and restitution 
coefficient 



Figures [3] (a) and (b) display, respectively, the mag- 
nified sequential plots of colliding clusters for a purely 
repulsive collision and a cohesive collision when the ini- 
tial temperature and the impact speed are T = 0.02e and 
V = O.Zy/eJm. From Fig. H we confirm that the con- 
tact duration for the cohesive collision is longer than that 
of the repulsive collision. [33j During the restitution, we 
also observe the elongation of the clusters along the z- 
axis in cohesive collisions, while we can not observe such 
a phenomenon in repulsive collisions. In both cases, the 
rotation of clusters is slightly excited after a collision. 

We firstly investigate the relation between the colliding 
speed and the restitution coefficient for a weak attractive 
case (c = 0.2). The cross points in Fig. [4] show the rela- 
tionship between the relative speed of impact scaled by 
the unit \J e/m, V = V/ \J e/m, and the restitution coeffi- 
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FIG. 3: (Color online) Sequential plot of collisions for (a) 
c=0.0 and (b) c=0.2 at i = 80, 85, 90, and 95, where t is time 
scaled by unit a / \fejm. 
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FIG. 4: Relationship between colliding speed and restitution 
coefficient for c = 0.2. V is relative colliding speed scaled by 
unit sJTJm. 



cient e. Sample average is taken over 100 different initial 
conditions for each speed and the error bar shows the 
standard deviation. From Fig. 0] we find that the resti- 
tution coefficient has a peak around the colliding speed 
V = 0.2y/e/m, which is attributed to the reduction of e 
caused by the attractive force in the lower impact speed. 
This tendency can also be observed in cohesive collisions 
of macroscopic bodies. 0, [24[ 

The solid line in Fig. 0] represents the theoretical 
prediction of cohesive collisions between viscoelastic 
spheres. [23j Here, we briefly summarize the theory of co- 
hesive collisions in Ref. 231 ]. Let us consider a head-on 
collision between elastic spheres of radii Ri and R2, each 
of which has mass of Mi and M2, respectively. The basic 
idea of their theory is to solve the time evolution equation 



of the deformation £(t) of the colliding spheres: 
tf(t) + F(£(t)) = 0, 

e(o) = o, m = v, 

where fj, is the reduced mass /i = (1/Mj + I/M2) -1 . 
is described as the function of the radius of contact 



a as 



£(«) = 



wxth R eff = f — 



R eff 

so that Eqs.((4|) are rewritten as 

£"(«)■ 2 
fia + /i ; a 



1 
~R~2 



0, 



«(0) = a imt , 6(0) =V -p| 



(4) 



area 



(5) 

(6) 
(7) 



where the prime denotes the differentiation with respect 
to a. We adopt a,i n it — (87rZ?7i?^^/3) 1 / 3 which is the 
contact radius of the bottom plane of the upper cluster 
with 7 ~ 0.104e/cr 2 estimated from the calculation of 
the attractive interaction between two clusters (see Ap- 
pendix ). We also estimate D as D = 3.28 x 10~ 3 a 3 /e 
from Young's modulus Y — 4:5<ie<T~ 3 and Poisson's ra- 
tio v — 7.74 x 10~ 2 which are obtained from another 
simulation. [361 ] 

They assume that the force F(a) between cohesive 
spheres comprises three kinds of forces: elastic force 
Fh(cl) characterized by Hertzian contact theory [I7I ITi!]. 
dissipative force Fdis (a) [6( , and cohesive Boussinesq force 
Fb{o) derived from JKR theory. [4oj] Thus, the total force 
can be expressed by 



F(a) = F H {a) - F B {a) + F dis (a). 



(8) 



Here, the sum of the elastic force and the Boussinesq 
force is given by 



F H {a)-F B (a) 



R, 



,3/2 



(9) 



with the surface tension 7 and D = (3/4){(l — v1)/2Y\ + 
(1 — ^|)/2l2} with Poisson's ratio Vi and Young's mod- 
ulus Yi for the cluster i — 1,2. Following the idea in 
Ref. [23| , we assume that the dissipative force is given by 

F dls (a) = Aa^-(F H (a) - F B (a)), (10) 
oa 

where A is a fitting parameter. 

We solve Eq. ([6]) with the initial speed V ranging from 
O.Ol^/m) 1 / 2 to 0.6(e/m) 1 / 2 by the use of the fourth order 
Runge-Kutta method to obtain the rebound speed which 
is the speed when the contact radius a becomes less than 

a sep = (^> j nD"fR 2 e ^/2 S j .[23| From the rebound speed 
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for each impact speed, we obtain the relationship be- 
tween the restitution coefficient and the impact speed. 
In Fig. [U we use A — O.lcr-v/m/e to draw the theoret- 
ical curve. In V < 0.2^/e/Vn, the discrepancy between 
our numerical results (cross points) and the theoretical 
result is large, which may be attributed to the rotational 
rebounds of clusters after collisions. On the other hand, 
the theoretical curve reproduces the results of simulation 
which excludes rotation of clusters (open circles) as will 
be explained later. 

Here, let us briefly comment on the dependence of 
the relative angle 6 Z on the numerical results. We have 
checked 9 Z dependence of the restitution for purely re- 
pulsive collisions, i.e. c — 0.0 at T = 0.02e. Figured] 
shows the relationship between impact speeds and resti- 
tution coefficients for 9 Z = 0, n/6, tt/3, and it/2. This fig- 
ure indicates that the relation between the impact speed 
and the restitution coefficient is not largely affected by 
the initial orientation, although the orientation around 
other axes may affect the relation. Thus, we will analyse 
only the results obtained with the fixed initial orientation 
9 Z = 0. 
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FIG. 5: Relation between relative colliding speed and resti- 
tution coefficient for different initial angles for c = 0.0 and 
T = 0.02e. 



Distribution of restitution coefficient 

In this subsection, we investigate the frequency dis- 
tributions of restitution coefficients for purely repulsive 
collisions and cohesive collisions, respectively. Figure 
[fDJa) shows the histogram of the restitution coefficient 
for purely repulsive collisions (c = 0.0). To obtain this 
result, we take 5000 samples at the fixed impact speed 
V = 0.02^/e/m. From Fig. [DJa), the frequency distri- 
bution can be approximated by the Gaussian (solid line) 
for purely repulsive collisions. 

On the other hand, Fig. [6^b) shows the frequency dis- 
tribution of the restitution coefficient for cohesive colli- 
sions (c = 0.2). To obtain this result, we take 995 samples 
at the fixed impact speed V = Q.l\Je/m. In Fig. HJb), 




0.1 0.3 0.5 0.7 0.9 1.1 

e 

(b) 

FIG. 6: Histograms of restitution coefficients for (a) c = 0.0, 
V = 0.02y/e/m, and (b) c = 0.2, V = O.ly/e/m. The solid 
line in (a) is the Gaussian fitting of the data. 



we find the existence of the two peaks around e = 0.448 
and e = 0.656, respectively, except for the main peak 
around e = 0.982. From the check of simulation movies, 
the collisions around these small peaks are produced by 
rotations after the collisions, while the most of bounces 
are not associated with rotations in the vicinity of the 
main peak around e = 0.982. It is reasonable that the 
excitation of macroscopic rotation lowers the restitution 
coefficient. 

Here, let us make another comparison of our simulation 
result with the theoretical curve drawn in FigJH The 
open circles in Figf?] are the mean values obtained by 
the data around the main peak for each impact speed to 
remove the effects of rotational bounces. It is obvious 
that the theory has a better fitting curve of the data 
when we remove rotational bounces. 

We shall comment on the fitting function of the main 
peak. FigureEJb) shows that the distribution around the 
main peak has an asymmetric profile, so that the Gaus- 
sian function may fail to fit the tail parts of the main 
peak. Figure [7] shows the semi- log plot of the simulation 
data around the main peak, where F(e) is the frequency 
of e. The reasonable fitting curves are represented by the 
solid lines, where lnF(e) = (19.6 ± 1.9)e+ (-14.7 ± 1.8) 
for e < 0.982 and lnF(e) = (-127 ± 28)e + (130± 28) for 
e > 0.982, respectively. Thus, the distribution of restitu- 
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FIG. 7: Semi-log plot of the main peak in Fig[SJb). The solid 
lines are the double exponential functions to fit the data. 



tion coefficients can be approximated by a combination 
of the double-exponential functions when the rotation is 
not excited after impacts. A similar tendency has also 
been observed in a recent experiment and a recent simu- 
lation of macroscopic collisions. [4~fl | 

Phase diagram of restitution coefficient 



tion of the impact speed under the fixed cohesive parame- 
ter c = 0.2 (Figfg](b)). Here, we find that the probability 
to emerge the modes (i) and (ii) decreases with the in- 
crease of the impact speed. In addition, the anomalous 
impact can be observed within the range of impact speed 
0.02 < V/ie/m) 1 ' 2 < 0.1. It is interesting that Figl^b) 
for V < 0.0Ay/e/m is almost the mirror symmetric one 
of FiglH^a) for c < 0.2, which suggests that the cohesive 
parameter plays a role of the impact speed. 
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FIG. 8: Probability diagrams classified by four collision 
modes with fixing (a) V = 0.02-^/e/m and (b) c = 0.2, re- 
spectively. 



As discussed in the previous subsection, some samples 
of the restitution coefficient exceeds unity even for co- 
hesive collisions. We can guess that most of colliding 
clusters coalesce when we use the collisional model with 
c = 1. Thus, it is important to know what process actu- 
ally occurs after a collision when the impact speed or the 
cohesive parameter c is given. In this subsection, we in- 
vestigate the emergence probability of four modes of the 
collisions (i) coalescence, (ii) bouncing, (iii) normal colli- 
sion with e < 1, and (iv) anomalous collision with e > 1. 
The coalescence (i) and the bouncing (ii) can take place 
only when the attractive interaction between the collid- 
ing clusters exists. Indeed, the bouncing occurs as the 
result of trapping by the potential well, if the rebound 
speed is not large enough. [42j 

Figure[5] (a) shows the phase diagram which is obtained 
under the fixed colliding speed V = 0.02^/e/m, where P 
represents the probability to observe each mode. This 
phase diagram shows that the regions for the modes (iii) 
and (iv) decrease with the increase of c. In the strong at- 
tractive case, c > 0.6, we cannot observe rebound modes 
(ii), (iii), and (iv). Here, we find that the anomalous im- 
pact can be observed for cohesive collisions with c < 0.4. 

We also categorize collisions into four modes as a func- 



Here let us reproduce the results of our simulation 
qualitatively by a phenomenology. Purely repulsive colli- 
sions, as we expect from Fig.JB^a), the probability density 
distribution of rebound speed V can be approximated by 
a Gaussian function 



p(V) = 



1 



'ira 



: exp 



(V - V m f 



(11) 



Thus, for given impact speed V, we can use Eq. fTTj) 
where V m and a in Eq.( are fitting parameters. 
Then, we calculate the probability to exceed the es- 
cape speed V* [43| from the attractive potential field 
$(r) = — 4ec(r / <t)~ 6 under the given V' . Here the es- 
cape speed V* of a rebounded cluster may be given by 



V*/ ^JTJm = ^2^(x*/a)/a^, 



(12) 



where x* — (2/c) 1 ^a at which the potential takes the 
minimum value. For example, the escape speed becomes 
V* = 0.015(e/m) 1 / 2 in the case of c = 0.2. 

Thus, from integrating the probability densities of V , 
the probabilities to observe modes (i), (iii), and (iv) are 
respectively given by 
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p(i) = 
p(iii) _ 
p(iv) 



pi V )dV 



-oo 
V 



p(V )dV 



2 1 I */5 



P (0, 



(13) 
(14) 
(15) 




0.2 0.4 0.6 0.8 1.0 

c 

(a) 



0.03 0.07 0.11 

V 
(b) 



FIG. 9: Probability diagrams from our theoretical argument 
for (a) V = 0.02y^7m and (b) c = 0.2. 



where erf (a;) is the error function erf(x) = 
J2 exp(—t 2 )dt. Here we ignore the distinction 
between the mode (i) and the mode (ii), because the 
most of bouncing clusters eventually coalesce after some 
numbers of collisions. 

Figures[9] (a) and (b) show the probability diagrams ob- 
tained from Eqs. gg]), gj}, and |gli)|. To draw FigG^b), 
we adopt V* = 0.018^/e/m which is slightly larger than 
the calculated value V* = 0MhyJe/m by Eq.(fl2]) for 
c = 0.2. Our phenomenology qualitatively reproduces 
the diagrams obtained by the simulation as in Figs [5] 
(a) and (b), although there are some quantitative dif- 
ferences between the simulation and the phenomenology. 
Indeed, the probability to appear the mode (iv) in the 
phenomenology decreases with the increase of Vyejm 
as in FiglQ] (b) , but this tendency cannot be observed in 
the simulation in Fig|8] (b) . 



DISCUSSION 

Let us discuss our results. We, mainly, discuss how 
the restitution coefficient depends on the cluster size in 
this section. Figure [TU] shows the relationship between 
the relative colliding speed of clusters and the restitution 
coefficient with different sizes of 236 atoms (G236), 433 
atoms (C433), and 682 atoms (C682), respectively, where 
we use the data obtained by the fixed parameters c = 0.2 
and T = 0.02e. As can be seen in Fig. [TU1 the restitution 
coefficients satisfies the scaling in which e(R/a) 317 is a 
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FIG. 10: Relation between colliding speed and restitution 
coefficient for clusters C236, C433, and C682- The restitution 
coefficients are scaled by R°' i17 , where R is the radius of each 
cluster. 



universal function of the impact speed, where R is the 
radius of each cluster. To obtain the scaling exponent, 
we first calculate the standard deviation for each rebound 
speed under the fixed value of the exponent. Next, we 
search the value of the exponent such that the maximum 
value of the standard deviations has a minimum value. 
To draw the solid curve in Fig[TUl we solve eq.© with 
the fitting parameter A = O.la^m/e for C682 with the 
aid of its radius R e ff — 3.23o\ This is interesting finding 
from our simulation which has not been predicted by the 
quasistatic theory of cohesive collisions for macroscopic 
bodies. 

From a simple phenomenology, we can understand that 
the restitution coefficient depends on the radius. How- 
ever, the phenomenology predicts that e(Rfa) 1 ^ 2 satisfies 
a scaling relation (see Appendix ). The discrepancy be- 
tween the phenomenology and the numerical observation 
indicates that our over-simplified theory is insufficient. 
We will need a more sophisticated theory to explain the 
exponent. 

We also simulate collisions between larger clusters than 
C682 by the use of C1055 and Ci466- Figure [TT1 (a) is the 
relationship between the impact speed and the restitution 
coefficients in the case of c = 0.0 under the initial tem- 
perature T = 0.02e. The squares, plus points, and circles 
show the averaged data of C682, C1055, C1466, respec- 
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tively. We take 10 samples for both C1055 and C1466 while 
100 samples for C682- Here we do not find any systematic 
relationship between the impact speed and the restitution 
coefficients in the cases of C1055 and Ci466- This can be 
attributed to the surface instability of the clusters arising 
from the weak attraction between "atoms" . 
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by the using the surface coated nanoclusters. For exam- 
ple, it has been demonstrated that hydrogen coated Si 
nanoparticles exhibit the weak attraction by H atoms on 



the surface. [34J, 135| We believe that our model captures 
the essence of such a system. For realistic simulations, 
we may need to carry out another simulation of the col- 
lision of H-passivated Si clusters by introducing suitable 
empirical potentials. As an additional remark, we should 
indicate that it is difficult to control the colliding speed 
and the initial rotation of the cluster in actual situations 
because the macroscopic motion of one cluster is also af- 
fected by thermal fluctuations. 

CONCLUSION 



FIG. 11: (a) Relation between impact speed and restitution 
coefficient in cases of C682, C1055, and Ci466- (b) Time evolu- 
tion of internal temperature of cluster C1055 in its free flight 
after initial equilibration to T = 0.02e. 



It is known that the instability of the spherical shape 
and the plastic deformation in a cluster cause the increase 
of the internal temperature of the cluster. [33l l34j We nu- 
merically performed free flights of cluster by the use of 
C1055 to check the time evolution of the internal tempera- 
ture of the cluster. FigureUJJb) shows the time evolution 
of the temperature inside the cluster C1055 after giving 
the translational speed V — 0.07y // 7Jm and the initial 
temperature T = 0.02e. Here we find the temperature 
increase during the free flight up to around T = 0.08e. 
Thus, we conclude that the maximum number of "atom" 
to reproduce the theory of cohesive collision is 682 in our 
system. 

We try to estimate the critical radius theoretically on 
the basis of the argument of capillary instability of elas- 
tic droplets (see Appendix ),[4J| but our over-simplified 
theory predicts that any elastic surface of spheres are 
unstable under the gravity. We should note that this cal- 
culation is based on theory of elasticity with zero shear 
modulus (Poisson's ratio is equal to -1). The calculation 
suggests that (i) we may not use theory of elasticity or 
(ii) zero shear modulus is unrealistic. We will, at least, 
need to discuss the capillary instability under the full set 
of elastic equations. From these arguments, we regard 
C682 as the maximum size to reproduce the quasi-static 
theory of cohesive collisions in our modelling. 

Although our simulation mimics impact phenomena of 
small systems subject to large thermal fluctuations, we 
should address that our model with small c may not be 
adequate for the description of many of realistic collisions 
of nanoclusters, where the cohesive interaction between 
clusters often prohibits the rebound in the low-speed im- 
pact. Namely, the corresponding value of the cohesive 
parameter is large in many of actual situations. How- 
ever, nanoscale impacts can be realized experimentally 



In conclusion, we have performed molecular dynam- 
ics simulations to investigate the behaviors of colliding 
clusters and the relationship between the restitution co- 
efficient and the impact speed. The results of our simula- 
tions have revealed that some aspects of the relationship 
can be understood by the quasi-static theory for cohesive 
collisions. [23J In addition, we have drawn the phase dia- 
gram of the restitution coefficient in terms of the impact 
speed and the cohesive parameter and explained them 
by a simple phenomenology. To clarify the mechanism 
of the emergence of the anomalous impact, it may need 
further investigation about the internal state of clusters 
during collision such as stress and modal analyses. 
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CALCULATION OF SURFACE TENSION 

In this appendix we explain how we calculate the sur- 
face tension 7 used to draw the theoretical curve in 
Fig. |H [2l| Let us assume that two identical clusters are 
in plane-to-plane contact with each other. When those 
clusters are located by the separation d, the surface en- 
ergy per unit area W is given by 21 1 



W 



B 



Undl 



,4 

d 2 



(16) 
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where B is the Hamaker constant B = A-K 2 €ca & p 2 with 
the cohesive parameter c and the number density p of 
each cluster. In our model, we use the number density 
becomes p = 0.4ct~ 3 , and d — 0.4o\ 

The surface tension 7 is equal to the energy per unit 
area to separate the two contacting plane to infinity. 
Thus, we obtain 7 as 

B 



7 



247rd2 



0.0261 x % = 0.1044-^. (17) 



DEPENDENCE OF e ON R 

Let us derive a scaling relation between the restitu- 
tion coefficient and the radius of cluster. Our assump- 
tion is that (i) the energy dissipation during a collision is 
originated from the sum of of the viscous force and the 
Boussinesq force, (ii) energy dissipation from the Boussi- 
nesq force is approximately given by the work during the 
detachment process of two coalesced clusters. 

Let us first estimate the energy dissipation caused by 
the surface tension. A pair of colliding clusters is par- 
tially coalesced as shown in Figure [HI where we assume 
that the deformation of two spheres are negligible, and 
contacted state can be characterized by a simple cut of 
the deformed region. Let 9 be the angle around the cen- 
ter of the upper sphere ranging from — 9q to 9$ under 
the assumption a small 9q. From this figure, the sur- 



where r is time scale of the dissipation. From the combi- 
nation of two terms in eqs. (fl~8)) and (p~9|) . we obtain the 
expression of the total energy loss during a collision 

Ei oss = E%: s ~W 7 ~ i? 2 ( Cl l/ n / 5 - c 2 V 4 / 5 ), (20) 

where c x ~ P 3 /5 tY 2 ^ and c 2 ~ (p /Y) 2 / 5 . Since Eq.® 
should be balanced with p R 3 V 2 (l — e 2 ), we obtain 

i?(l-e 2 ) ~ — {dV 1 ^ 5 - c 2 V- 6 ^). (21) 
Po 

Thus, our phenomenology suggests that R}l 2 e is inde- 
pendent of the radius of the colliding spheres. 



INSTABILITY OF AN ELASTIC DROPLET 

In this appendix we investigate the instability of the 
surface profile of clusters on the assumption that the in- 
ternal modes of the cluster are expressed by those of an 
isotropic elastic sphere. When the shear stress can be 
ignored, the stress tensor aij can be written as 



KV ■ uSi 



(22) 



where K is the bulk modulus, <5y is Kronecker delta, and 
u is the strain. Thus, the equation of motion in the bulk 
becomes 



pu = V(KV ■ u) = -Vp, 



(23) 





e\ ( 


J. /> 









FIG. 12: Schematic figure of contacting identical spheres. 



where p is the density and p = —KV ■ u is the effec- 
tive pressure. Thus, the problem can be mapped onto 
a problem of perfect fluid. Thus, the dispersion relation 
is linearized equation R(9,<fi,t) = R a + £(6,<fr,t) can be 
written as 



7 ?(Z-l)(Z + 2) 
P Rl 



(24) 



as in the case of a liquid droplet [44( , where I is the index 
of Legendre polynomial. 

When we introduce the gravity in this perfect fluid 
model, the scalar potential <I> defined by v = V<i> satisfies 



face area of a cut hemi-sphere is approximately given by 
ttR 2 9q — 7ri?£ with the deformation £ ~ 9^R. Thus, the 
work needed to pull off two spheres is 

2/5 



W 1 ~ 27T7.R&, 



R 2y4/5 



(18) 



where we use the estimation £, m ax oc 
(p Q R 3 /Yy/R) 2 / 5 V' l/5 - (po/Y) 2 / 5 RV i/r ° on the basis of 
the theory of elasticity, where po is the density. [18 1 

On the other hand, the energy dissipation of repulsive 
spheres is given by@ 



Ef:: s ^pi /5 TY 2 ^R 2 v n / 5 , 



(19) 



d t <f> + P+-v 2 +gz = f(t), 



(25) 



where where dt is the time derivative, P — Jdp/p(p), 
and f{t) is an arbitrary function of time, g and z are the 
gravitational acceleration and the relative vertical posi- 
tion from the center of mass of the sphere. Choosing f(t) 

satisfying f(t) = p + 7 f-^- + -^j with the surface ten- 
sion 7, curvatures R\ and R 2 , Eq. ([23]) can be rewritten 
as 
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where A(M) = - {^M ( sil ^Jj) + 
derive Eq. (f2"6| we have used $ = $, 
and £ = v r = d<fr/dr at r = Rq. 



l a- 1 

sin^ 9 d(j} 2 



u = V*, 



By using the expansion $ in terms of r l and the spher- 
ical harmonic function Yi tTn (9, (/>), we may obtain 



= _^ i(i _ 1)(i+2) _|. 



(l-m)(l + m) (l-m+l)(l + m+l) 



(21- 1)(2/ + 1) 



(2Z + 1)(2Z + 3) 



(27) 



r 



where we have used the formula 



cos 9Y lm = \ \ fni i i oX Yi+i,m 



(2/ + l)(2/ + 3) 



(I — m)(l + m) 
(21- 1)(2/ + 1) 



Yl-l.m- (28) 



Therefore, u> n ,i,m becomes complex, if the radius exceeds 
the critical radius 



R 



(l,m) 
0,c 




(/ — m)(Z + m) 



1)(2Z 




(Z-m + l)(Z + m + l) 
(2Z + 1)(2Z + 3) 



1/2 



(29) 



Equation (|29|) implies that the mode with I = 1 is always 
unstable for the perturbation. Thus, we conclude that an 
accelerated elastic sphere is unstable, which is similar to 
the instability of a raindrop of the perfect fluid because 
the neutral mode I = 1 does not have any recovering 
force. 
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